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ABSTRACT 

Large telescopes have allowed astronomers to observe galaxies that formed as early as 
850 million years after the Big Bang. We predict when the first star that astronomers 
can observe (i.e., in our past light cone) formed in the universe, accounting for the first 
time for the size of the universe and for three essential ingredients: the light travel time 
from distant galaxies, Poisson and density fluctuations on all scales, and the effect of 
very early cosmic history on galaxy formation. We find that the first observable star is 
most likely to have formed 30 million years after the Big Bang (at redshift 65). Also, 
the first galaxy as massive as our own Milky Way likely formed when the universe 
was only 400 Myr old (at redshift 11). We also show that significant modifications 
are required in current methods of numerically simulating the formation of galaxies 
at redshift 20 and above. 
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1 INTRODUCTION 

The formation of the first stars, in a universe that had 
not yet suffered chemical and magnetic stellar f eedback, 
is a well-defined problem for theorists (Barkana. fc Loebl 
I2OOII ). The advent of large telescopes has accelerated the 
observational quest for the first objects as astronomers have 
detected the ligh t emitted by increasingly distant galaxies 
(|Fan. et al.l 120031 ). A theoretical understanding of the first 
stars is of great importance as observations approach the 
pristine regime of the early universe. 

Stars are understood to form at high redshift out of gas 
that cooled and subsequently condensed to high densities in 
the cores of dark matter halos. Since metals are absent in the 
pre-stellar universe, the earliest available coolant is molec- 
ular hydrogen {H2), and thus the minimum halo mass that 
can form a star is set by requiring the infalling gas to reach 
a temperature > 1000 K requir ed for exciting the rot ational 
and vibrational states of H2 IjTegmark et al.l 119971 ) . This 
has been confirmed with high-resolution numerical simula- 
tions containing gravity, hyd rodynamics, and chemical pro- 
cesses in the primo rdial gas (Ab el. Bryan fc Normarj|2002l : 
iBromm. Coppie fc Larson 2002 ). These simulations showed 
that the first star formed within a halo containing ~ 10^ Mq 
in total mass; however, as we show, they could not estimate 
when the first star formed. 
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2 LIMITATIONS OF SIMULATIONS 

In the real universe, the cosmological redshift z of each 
source implies a corresponding cosmic age at which the 
source is seen. Thus, the observable universe can be divided 
into thin spherical shells around us, each containing a por- 
tion of the universe that is observed at a given redshift z and 
is seen at the corresponding age. Galaxy formation in differ- 
ent shells is correlated due to density fluctuations on various 
scales. A simulation that wishes to determine when various 
objects are observed to form must encompass a sufficient 
portion of our past light cone. 

Even with adaptive grid codes, full simulations are lim- 
ited by current technology t o tiny regions and only form a 
first star below redshift 20 l|Abel. Bryan fc NormanI |2002|) 
or 33 for the most recent simulations ( OShea fc NormanI 
[2OO6). Even when great sacrifices are made in the real- 
ism, so that hydrodynamics and chemical evolution are 
neglected and the collapse of a star is only partially re- 
solved, simulations still cannot capture the proper cosmo- 
logical context. A direct simulation of the universe out to 
the spherical shell at redshift 70 would require a simulated 
box of length 25, 000 Mpc on a side [measured in comov- 
ing distance, which equals physical distance times a factor 
of (1 + z)]. The need to resolve each W^Mq halo into 500 
particles (requ ired to de termine the halo's evolution even 
crudely (Sprin gel fc Hern quist 200Ji)) determines the mass 
of each simulated particle, while the mean cosmic density 
of 4 X 10^°Mq/Mpc^ yields the total mass of the simu- 
lated box, and implies a required total of 3 x 10^^ parti- 
cles. The maximum number achieved thus far is only ^ 10^" 
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IjSpringel et alJbOOSl ). Early times have been pr obed by res- 
imulating regions from a large initial simulation (|Reed et alj 
12003 '). forming star- forming halos at 2 = 47 (the earliest 
to date), but with cosmological parameters for which the 
first galaxy actually formed at 2 ~ 82 according to our 
model below. Thus, the actual starting redshift of the era 
of galaxy formation can curre ntly only be pred i cted a nalyt- 
ically. For similar parameters, iMiralda-Escudg (120031) made 
a rou gh analytical estimate (using the IPress fc Schechteij 
|[197J) halo mass function) and foun d that the fir s t star 
formed at 2 



48. A rough estimat e by I Wise fc Abell ( 
found 2 ~ 71 using the improved ISheth fc TormerJ ( 



2005 ) 



19991 ) 



mass function (see also § 3). 

Simulations of the first stars also face a difficulty in 
the way their initial conditions are determined. Recent 
observations of the cosmic microwave background (CMB) 
l|Spergel et al.l I2OO6I ) have confirmed the notion that the 
present large-scale structure in the universe as well as the 
galaxies within it originated from small-amplitude density 
fiuctuations at early cosmic times. Simulations thus assume 
Gaussian random initial fiuctuations as might be generated 
by a period of cosmic infiation in the early universe. The evo- 
lution of these fluctuations can be calculated exactly as long 
as they are small, with the linearized Einstein-Boltzmann 
equations. The need to begin when fluctuations are still lin- 
ear forces numerical simulations of the flrst star-forming ha- 
los to start at a h igh redshift, the highest to-date being 600 
l|Reed et al.ll2005h . Current simulations neglect the contri- 
bution of the radiation to the expansion of the universe and 
assume the Newtonian limit of general relativity which is 
valid in simulation boxes that are small with respect to the 
horizon. However, for a halo that forms around redshift 65, 
we find below that the corresponding region already must 
have had a moderately non-linear overdensity of 5 ~ 20% at 
2 = 600 {8 is the density perturbation divided by the cos- 
mic mean density). Indeed, we find that for initial conditions 
that ensure 1% accuracy in 5, simulations must start deep 
in the radiation dominated era, at 2 > 10® . Furthermore, we 
show that a small error in 5 leads to a much larger error in 
the abundance of halos at each redshift. The non-linearity 
of the initial conditions represents a previously-unrecognized 
fundamental limitation of current methods of simulating the 
formation of the first star-forming halos. 



3 SPHERICAL COLLAPSE 

While numerical simulations cannot determine the forma- 
tion time of the first observable star, galaxy and halo forma- 
tion can also be understood and described with a standard 
analytical model that has been quantitatively checked with 
simulations where the latter are reliable. The standard ana- 
lytical model for the abundance of halos ([ Press & Schcchtor 
I1974I : iBond. Cole. Efstathiou fc Kaiseij|l99ll ) considers the 
linear density fiuctuations at some early, initial time, and 
attempts to predict the number of halos that will form at 
some later time corresponding to a redshift 2. First, the 
fiuctuations are evolved to the redshift 2 using the equa- 
tions of linear fiuctuations, i.e., they are artificially linearly 
extrapolated even if they may become non-linear in some 
regions. The average density is then computed in spheres 
of various sizes. Whenever the overdensity in a sphere con- 



taining mass M rises above a critical threshold 5c{M,z), 
the corresponding region is assumed to have collapsed by 
redshift 2, forming a halo of at least mass M. Thus, the 
model separates the linear grow th of fluctuations , which we 
calculate with standa r d codes (ISeliak fc Zaldarr iaga 1996|; 
iNaoz fc Barkanall2005l : lYamamoto. Sugivama fc SatOiil998l ) , 
from the non-linear collapse of objects, which determines the 
effective linear threshold 5c{M,z). 

Within this model, the non-linear collapse of a halo 
is analyzed by considering a uniform, spherically symmet- 
ric fluctuation whose collapse can be calculated by numer- 
ically solving ordinary differential equations; 5c{M,z) is 
then determined as the value of the linearly-extrapolated 
5 of the spherical region at the moment when the ac- 
tual, non-linear 5 diverges (i.e., the region collapses to a 
point). The classical calculation of spherical collapse in 
an Einstein de-Sitter (EdS) universe (i.e., a flat matter- 
dominated universe) yields 5c ~ 1.686, a value that is 
independent o f the halo mass and the collapse redshift 
IjCunn fc GottI Il972l ). Given the power spectrum of the 
initial fluctuations, together with their linear growth, we 
can determine the chance that regions containing various 
initial masses M will reach the threshold 5c at some fi- 
nal redshift 2, and this forms the basis for estimating the 
abundance of halos of mass M that form at a specific 2 
(|Pr ess fc Schechteij [l97i : iBond. Cole. Efstathiou fc Kaiser! 
il991i). We use th e rece ntly-modified form of this model 



Sheth fc TormenI (119991). with the u pdated parameters 
suggested by ISheth fc TormenI (|2002| )]. which fits halo 
abundances in numerical simulations very accurately in 
regi mes that include M = 10^^-10^^ Mg halo s at 2; = 0- 
10 (|jenkins et al.l I2OOII : ISpringel et all l2005l) and M = 
lO'^ Mp) halos at 2 = 15-30 IJBarkana fc Loebl |2004| : 



lY oshida. et al.ll2003|). We note that other suggested mass 
functions (^e.g.. | Reed et al.ll2006l ) are in agreement with the 
ISheth fc TormenI 1 19991 ) mass function if the updated pa- 
rameters are used. Based on the plausible physical argu- 
ments behind this model together with its quantitative suc- 
cess out to redshift 30, we extrapolate it to predict halo 
abundances at still-higher redshifts. 

Our results assume cosmological parameters match- 
ing the latest CM B and weak lensing observations 
(jSpergel et al.M2006l ). For the contributions to the energy 
density, we assume ratios relative to the critical density of 
Q.m = 0.299, f^A = 0.701, and fib = 0.0478, for matter, 
cosmological constant, and baryons, respectively. We also 
assume a Hubble constant Hq = lOO/i km s~'^Mpc~^ with 
h = 0.687, and a primordial power-law power spectrum with 
spectral index n = 0.953 and as = 0.826, where ag, is the 
root-mean-square amplitude of mass fluctuations in spheres 
of radius 8 /i~^ Mpc. 

Within this analytical model, we can account for early 
cosmic history when predicting the abundance of high- 
redshift stars and galaxies. We flrst calculate the correct 
value of Sc{M,z) with a spherical collapse calculation that 
includes the full formation history of each halo (Figure 1), 
starting out with a small perturbation and following its grav- 
itational evolution until it collapses. We begin our calcu- 
lation at a very high redshift when 5 <^ 1, which implies 
that the overdense region containing the mass M is larger 
than the horizon. In this regime we use the cosmological 
Friedmann equations of general relativity, applied to the 
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Figure 1. Evolution of the fractional overdensity 8 for a spherical 
region containing 10^ Mq that collapses at 2 = 66. We show the 
fully non-linear & (solid curve) and the linearly-extrapolated <5 
(dashed curve). We indicate the redshifts of halo collapse (zcoll)i 
cosmic recombination (zroc): matter-radiation equality (zeq), and 
entry into the horizon (zcntcr). 
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Figure 2. Cumulative halo mass function n(> M). At 2 = 66 
(top panel) we compare our full result for n{> M) (solid curve) to 
that resulting from adopting the EdS value of 5c = 1.686 (dashed 
curve). At 2 = 20 (bottom panel), we show the relative error in 
n(> M) that results from using 5c = 1.686 instead of our full 
result. 



overdense region that we are following. We begin in the 
radiation-dominated era, and include throughout the con- 
tribution of the radiation to the cosmic expansion. Once the 
fluctuation enters the horizon, the overdense dark matter 
continues to collapse due to gravity, but the radiation pres- 
sure suppresses the sub-horizon perturbation in the photon 
density, and the coupling of the baryons to the photons (via 
Compton scattering) keeps the baryon perturbation small as 
well. Thus, in this regime we neglect the perturbation in the 
radiation and baryons and continue to evolve the spherical 
collapse of the dark matter perturbation. The formation of 
hydrogen at cosmic recombination {z ~ 1200) decouples the 
cosmic gas from its mechanical drag on the CMB, and the 
baryons (which constitute a significant 17% of all the mat- 
ter) subsequently begin to fall into the pre-existing gravi- 
tational potential wells of the dark matter. In this regime 
we evolve the coupled collapse of the dark matter and the 
baryons, both evolving under their mutual gravity but with 
different initial conditions. We find that a halo destined to 
contain the first observable star in the universe reaches a 
5 ~ 1% at 2 ~ 10® and grows to 5 ~ 6% at matter-radiation 
equality and 5 ~ 13% at cosmic recombination. Since the 
fiuctuation enters the horizon very early in the radiation- 
dominated era, the final value of Sc{M, z) for a given M and 
z is insensitive to the precise starting redshift (as long as 
it is much higher than equality) or the precise treatment of 
the fiuctuation's horizon crossing. 

We find that 5c is essentially independent of M, and 
for our assumed parameters it is lower than the EdS value 
by ~ 1% X (1 -I- z)/2Q in the range z = 9-80. Adding the 
contribution of the radiation to the expansion of the uni- 
verse, as well as adding the baryon-photon coupling to the 



collapse, both reduce the extrapolated linear perturbation 
at collapse &c compared with the EdS value. These added 
physical effects reduce the collapse efficiency since only the 
dark matter component takes part in the collapse through- 
out. The radiation is kept smooth by its own pressure and 
it never participates in the collapse, while the baryons only 
gradually recover from their coupling to the photons and 
begin to catch up to the dark matter; during this recov- 
ery stage, as long as the baryon perturbations remain much 
smaller than those of the dark matter, the baryons essen- 
tially do not participate in the collapse. Now, any reduction 
in the matter fraction that collapses (compared to 100% 
in the usually- assumed case of pure cold dark matter) de- 
presses the linear evolution of the density perturbation more 
strongly, while the non-linear perturbation is larger and is 
thus less affected by the components that do not help in the 
collapse. Therefore the linear perturbation reaches a lower 
value of 5c when the non-linear perturbation collapses. Note 
that while the change in 5c niay seem small, when dealing 
with rare halos, even a change of a few percent in 5c can 
change the halo abundance at a given redshift by an order 
of magnitude (Figure 2). Even at redshift 20, the 1% change 
in 5c compared to its previously-assumed value changes the 
cumulative halo mass function n{> M) by > 10%. 



4 STATISTICAL CONSIDERATIONS 

As noted earlier, the minimum halo mass that can 
form a star at high redshift is determined by H2 
cooling (ITegmark et al.l 19971 ). Numerical simulations 
IIAbel. Bryan fc NormanI I2OO2I: iBromm. Coppie fc LarsonI 
l2002l : iFuller fc CouchmanJ I2OO0I : lYoshida. et all l2003l : 
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iReed et al.ll2005l ) imply a minimum circular velocity Vc ^ 
4.5 km/s, where Vc = \JGM/R in terms of the halo virial 
radius R. This yields a corresponding minimum mass at each 
redshift z' , Mmin(z')- The simulations show that in a halo 
above the minimum mass, the gas cools in the dense center 
and forms at least one star very quickly; this is understood 
theoretically since both the cooling time and the dynamical 
time are a small fraction of the cosmic age at that time. 

Given the comoving number density of such star- 
forming halos at each redshift, n{Mmi-a{z'),z'), the mean 
expected number of star-forming halos observed at redshift 
z or higher is 



{N){z) = 



n{Mn 



(z),z) — dz , 



(1) 



where dV is the comoving volume of the spherical shell in 
the redshift range z' to z' + dz' . Our results are independent 
of the adopted majdmum redshift as long as imax > 80. 
Although the mean expected number {N){z) is determined 
theoretically by the cosmological parameters, in practice our 
observable universe presents us with a single instance and 
the actual observed redshift is subject to Poisson and density 
fluctuations. 

The effect of Poisson fluctuations is easily calculated. 
The chance of observing at least one star-forming halo 
above redshift z with Poisson fluctuations included is 1 — 
exp[— (A'^)(2:)], which yields a spread in the redshift of the 
first such halo of Az ;« 3 at 1 — cr. In addition, the effect of 
density fluctuations on all scales is included statistically in 
our analytical model for the mean halo abundance at each 
redshift, but the correlations in the halo abundance between 
nearby points are not properly included. In other words, our 
calculation assumes that the relevant redshift shell gives us a 
fair sample of the density fluctuations on all scales. We first 
note that theoretically we expect that if we included the fully 
correlated density field, the effect on the predicted redshift 
of the first star-formation halo would be much smaller than 
the spread we find due to Poisson fiuctuations. Since the 
universe is homogeneous on large scales, large-scale modes 
can shift the halo abundance {N){z) only by a small Az; in- 
deed, onl y 10 Mpc scales and smaller can produce a Az > 1 
(jBarkana fc Loeb 2004 ). Such small-scale modes are indeed 
very well sampled within the spherical redshift shells that 
have a very large radius of ^ 12, 500 Mpc. We have also 
verified this with a direct quantitative test. We have used 
the power spectrum to produce instances of the full, cor- 
related density field in three-dimensional boxes. We have 
then modified the abundance n(Afniin(2), z), increasing it in 
high-density regions and decreasing it in low-density regions 
according^to a model that fits the results of numerical sim- 



ulations ( Barkana fc Loebll2004l ). In a 1 Gpc'^ box resolved 
into 256'^ cells, we found that including a correlated density 
field shifts the typical redshift of the first star-forming halo 
in the box by a Az < 1, which is already a smaller effect 
than the Poisson redshift spread. The spherical shells in the 
real universe have a larger volume by a factor of ~ 2000, so 
the density fluctuations are far better sampled than in our 1 
Gpc^ box, and the effect of the correlations on the redshift 
of the first star-forming halo in the universe is negligible. 
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Figure 3. The probability of finding the first star as a function 
of redshift. Vertical lines show the median redshift as well as the 
central 1 — a (68%) range and 2 — a (95%) range. 



5 PREDICTIONS 

The most likely redshift for the first observable star is 65.4 
(Figure 3), while the median redshift at which there is a 50% 
chance of seeing the first star is 65.8, which corresponds to 
a cosmic age of t = 31 Myr, less than a quarter of one per- 
cent of the current age of 13.7 Gyr. The 1 — a (68%) range 
is z = 64.7-67.3 (or t = 30-32 Myr) and the 2 - cr (95%) 
range is z — 63.9-69.4. Note that even an error in the min- 
imum Vc as large as a factor of 1.5 would cause only a 9% 
change in the expected formation redshift. Note also that 
in terms of the total cosmic mass fraction contained in such 
halos, the first star-forming halo corresponds to an 8.3 — a 
density fluctuation on the mass scale of 10^ M0. Future sim- 
ula tions will allow us to ch eck the validity of extrapolating 
the lSheth fc TormenI l|l999l ) halo mass function to such rare 
fluctuations. 

We can consider a similar question for other popula- 
tions of halos (Figure 4). For example, the radiation from 
the first stars is expected to eventually dissociate all the H2 
in the intergalactic medium, leading to the domination of 
a second generation of larger galactic halos where the gas 
cools via radiative transitions in a tomic hydrogen and he- 
lium (|Haiman. Rees. fc Loebl 119971 ). Atomic cooling occurs 
in halos with Vc > 16.5 km/s, in which the infalling gas is 
heated above 10,000 K and is ionized. We find that the flrst 
galaxy to form via atomic cooling formed at z = 46.6lQg, 
which corresponds to t = 52^2 Myr, where we have indicated 
the median values and 1 — cr ranges. We also predict that 
the first halo as large as that of our own Milky Way galaxy 



(|Battaglia et al.ll2005l ) formed at z = ll-lia2. or t = 408 
Myr. Al so, the first halo with t he mass of the Coma galaxy 
cluster l|Lokas fc Mamonll2003l ) formed at z = 1.24lo;Jo, or 
t = 5.0 ± 0.4 Gyr. 
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Figure 4. The median redsiiift for the first appearance of various 
populations of halos. We consider halos above a minimum circular 
velocity (left panel) or minimum mass (riglit panel). We indicate 
in particular the first star-forming halo in which H2 allows the 
gas to cool, the first galaxy that forms via atomic cooling (H), as 
well as the first galaxy as massive as our own Milky Way and the 
first cluster as massive as Coma. The horizontal lines indicate the 
elapsed time since the Big Ba ng. We compare the results with the 
cosmological parameters from lSpergel et al.l 1120061) [solid curve] to 
those with the parameters from lViel et al.l l l2006r [dashed curve]; 
the difference represents roughly the systematic la error due to 
current uncertainties in the values of the cosmological parameters. 



always smaller than the age of the universe at the time by 
at least an order of magnitude. Since the baryonic compo- 
nent collapses by another factor of ~ 20 until the collapse is 
halted by angular momentum (assuming a typical halo sp in 
parameter of ~ 5%; see, e.g.. IMo. Mao, fc Whi"t3 llQQSh ). 
the dynamical time within a disk or a star-forming region is 
likely to be shorter still by another two orders of magnitude. 
In summary, by accounting for the size of the observable 
universe as well as the early cosmic history of individual ha- 
los (Figure 1), our results extend the era of star formation 
much earlier than current simulations suggest. An analyti- 
cal model of galaxy formation that has been normalized to 
simulation results allows us to predict the earliest time that 
various halo populations can be observed (Figure 4). Since 
the statistical uncertainty in this time is dominated by Pois- 
son fluctuations, we can analytically predict its full proba- 
bility distribution (Figure 3). Our results also show that in 
order to numerically simulate galaxy formation accurately 
at z ^ 20, current simulation methods must be revised to 
include the effect on a forming halo of the early cosmic his- 
tory of radiation and baryons. Simulations that neglect this, 
even if they overcome the statistical challenges by covering 
the volume of the entire observable universe, will underesti- 
mate the number of halos at 2 = 20 (Figure 2) by 9% for all 
halos above the H2 cooling threshold (6 x 10^ Mq at z = 20), 
and by 15% for all halos above the atomic cooling threshold 
(3 X IO^Mq at z = 20). 
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6 DISCUSSION 

There remains a ~ 10% systematic uncertainty in the above 
1 + z values due to current uncertainties in the cosmolog- 
ical parameters. Varying the cosmological parameters has 
a negligible effect on the critical overdensity (Sc) but sig- 
nificantly changes the power spectrum and also the spatial 
scale corresponding to a given halo mass or circular veloc- 
ity. We illustrate the uncertainties in Figure 4 us i ng th e 
alternate cosmological parameters from IViel et al.l l|2006l ): 
n™ = 0.253, ^A = 0.747, Qb = 0.0425, h = 0.723, n = 0.957 
and as = 0.785. These parameters also represent typical 
1 — a errors, in terms o f the parameters uncertainties given 
bv lSpergel et al.l (|2006l ). We obtained for these alternate pa- 
rameters that the first star formed ai z = 60.0, which rep- 
resents a 9% reduction in 1 -f 2 compared to our standard 
parameters; also the value of 1 -I- 2 is lowered by 9% for the 
first atomic-cooling halo, 7% for the first Milky Way-mass 
halo, and 8% for the first Coma-mass cluster. 

On the other hand, astrophysical uncertainties about 
metallicity and feedback can only affect the precise prop- 
erties of a halo of a given mass, but not whether or not it 
forms. Also, in a halo with a short cooling time, processes 
within the halo should not induce a significant uncertainty 
in the formation time of the first stars within it, since a 
halo virializes at a density of ~ 180 times the cosmic mean 
density, and so the gravitational dynamical time within it is 
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